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Abstract 

Quantum fluctuations are believed to play an important role in the 
thermalization of classical fields in infiationary cosmology but their rele- 
vance for isotropization/thermalization of the classical fields produced in 
heavy ion collisions is not completely understood. We consider a scalar (^"^ 
toy model coupled to a strong external source, like in the Color Glass Con- 
densate description of the early time dynamics of ultrarelativistic heavy 
ion collisions. The leading order classical evolution of the scalar fields 
is significantly modified by the rapid growth of time-dependent quan- 
tum fluctuations, necessitating an all order resummation of such "secu- 
lar" terms. We show that the resummed expressions cause the system to 
evolve in accordance with ideal hydrodynamics. We comment briefly on 
the thermalization of our quantum system and the extension of our results 
to a gauge theory. 



1 Introduction 

A remarkable outcome of the heavy ion experiments at RHIC [1-4] is the large 
elliptic flow observed in the collisions. Phenomenological hydrodynamical mod- 
els that fit the RHIC data appear to require that the quark gluon matter has 
a very small value for the dimensionless ratio of the viscosity to the entropy 
density [5]. This ratio 77/s, a measure of the "perfect fluidity" of the system, is 
estimated to be < 5/47r [6], where rj/s = l/47r is a conjectured universal lower 
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bound [7,8]. Its was shown recently [9 11] that the degree of perfectness of the 
quark-gluon fluid produced at RHIC is sensitive to details of the initial spatial 
distribution of the produced matter at the onset of hydrodynamic flow. 

An important feature of the hydrodynamic models is that they require very 
early thermalization after the collision. Estimates for the thermalization time, 
which range from Treiax ~ 0.6 — 1 fm [12-14], are difficult to reconcile with a sim- 
ple picture of thermalization arising from the rapid scattering of quasi-particles 
at rates greater than the expansion rate of the fluid. The uncertainty principle 
tells us that for Treiax < 1 &n, modes with momenta ~ 200 MeV are not even 
on-shell, let alone amenable to being described as quasi-particles undergoing 
scattering. While a quasi-particle description is not essential to thermalization, 
it is the simplest one, and other realizations are more complicated. With regard 
to the issue of flow however, it is sufficient to note that one requires primar- 
ily that matter be isotropic and (nearly) conformal to obtain a closed form 
expression for the hydrodynamic equations [15]. 

How isotropization and (subsequently) thermalization is achieved in heavy 
ion collisions is an outstanding problem which requires that the problem be 
considered ab initio. What this means it that one needs to understand and 
compute the properties of the relevant degrees of freedom in the nuclear wave- 
functions and how these degrees of freedom decohere in a collision to produce 
quark-gluon matter. An ab initio approach to the problem can be formulated 
within the framework of the Color Glass Condensate (CGC) effective field the- 
ory, which describes the relevant degrees of freedom in the nuclei as dynamical 
classical fields coupled to static color sources [16-18] . The computational power 
of this approach is a consequence of the dynamical generation of a semi-hard 
scale, the saturation scale [19,20], which allows a weak coupling treatment of the 
relevant degrees of freedom [21-23] in the high energy nuclear wavefunctions. 

There has been significant progress recently in applying the CGC effective 
field theory to studying the early time behavior of the matter produced in heavy 
ion collisions. Inclusive quantities such as the pressure and the energy density in 
this matter (called the Glasma [24]) can be written as expressions that factorize 
the universal properties of the nuclear wavefunctions (measurable for instance in 
proton-nucleus or electron-nucleus collisions) from the detailed dynamics of the 
matter in collision [25 27]. Key to this approach are the quantum fluctuations 
around the classical fields in the wavefunctions and in the collision. Quantum 
fluctuations that are invariant under boosts can be isolated in universal func- 
tional that evolve with energy. There are however also quantum fluctuations 
that are not boost invariant which are generated during the collision. These 
quantum fluctuations can grow rapidly and therefore play a significant role in 
the subsequent temporal evolution of the Glasma. 

The problem of how to treat these so-called "secular divergences" of per- 
turbative series is very general and occurs in a wide variety of dynamical sys- 
tems [28]. In particular, the role of time dependent quantum fluctuations in 
heavy ion collisions bears a strong analogy to their role in the evolution of the 
early universe [29]. In the latter case, quantum fluctuations around a rapidly 
decaying classical fleld, the inflaton, are enhanced due to parametric resonance, 
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and it is conjectured that this dynamics termed "preheating" [30] may lead to 
turbulent thermalization [31] in the early universe. 

It is therefore very important to understand the precise role of these quantum 
fluctuations in heavy ion collisions to determine whether they play an analogous 
role to that in the early universe in the isotropization/thermalization of the sys- 
tem. Their computation in a gauge theory is quite involved so for simplicity, we 
shall in this paper first attack this problem in a scalar field theory. Like QCD, 
the coupling is dimensionless in this theory and the fields are self interacting. 
In addition, we choose initial conditions for our study that are similar to those 
in the CGC treatment of heavy ion collisions. It must be said at the outset 
that there are important differences between the two theories and there is no a 
priori guarantee that the lessons learnt in one case will translate automatically 
to the other. 

The CGC initial conditions, for weak couplings <? <C 1, specifically lead to a 
power counting scheme where the leading contribution to inclusive quantities is 
the classical contribution of order 0(l/g^). Quantum corrections begin at 0(1) 
and their contribution can be expressed as real-time partial differential equations 
for small fluctuations in the classical background, with purely retarded initial 
conditions. We will show that there arc modes of the small fluctuation fleld that 
grow very rapidly and can become as large as the classical fleld on time scales of 
interest in the problem. We observe that there are two sorts of rapidly growing 
modes of the fluctuation fleld. One are modes that enjoy parametric resonance 
and grow exponentially. These modes are however localized in a rather narrow 
resonance band. The zero mode and low lying modes grow linearly and can 
also influence the temporal evolution of the system. Both sorts of "secular" 
terms can be isolated and resummed to all orders in perturbation theory. The 
resulting expressions are stable and can be expressed as an ensemble average over 
a spectrum of quantum fluctuations convolved with the leading order inclusive 
quantity which, for a particular fluctuation fleld, is a functional of the classical 
field shifted by that quantum fluctuation. We note that a similar observation 
was made previously in the context of inflationary cosmology [32 34] . 

The fact that one can express resummed expressions for the pressure and 
energy density as ensemble averages over quantum fluctuations has profound 
consequences. Without resummation, the relation between the energy density 
and the pressure is not single valued. For the resummed expressions, while 
the relation between the pressure and energy density is not single valued at 
early times, it becomes so after a flnite evolution time. This development of an 
"equation of state" therefore allows one to write the conservation equation for 
the resummed energy momentum tensor T'"' as a closed form set of equations, 
which are the equations of ideal hydrodynamics. This of course suggests that 
the system behaves as a perfect fluid. If the considerations in our paper can 
be applied to a gauge theory, the result would have significant ramifications for 
the interpretation of the heavy ion experiments and the extraction of rj/s in 
hydrodynamical models. 

The evolution of the system towards the equation of state characteristic of 
hydrodynamic flow can be interpreted as arising from a phase decoherence of 
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the different classical trajectories of the energy momentum tensor for different 
initial conditions given by the ensemble of quantum fluctuations. For a scalar 0^ 
theory, the frequency of the periodic classical trajectories is proportional to the 
amplitude. Therefore, for different initial values of the amplitude, the different 
trajectories are phase shifted. The ensuing cancellations between trajectories 
results in the single valued relation between the pressure and the energy density. 
While it appears that decoherence can arise from the zero mode and near lying 
modes alone, the inclusion of the resonant band significantly alters the deco- 
herence of the system. Similar behavior has been seen in models of reheating 
after inflation [35 37]. In particular, one sees that quantum de-cohcrcncc of 
the inflaton field leads to a transition from a dust-like equation of state to a 
radiation dominated era. 

It is interesting to ask whether the decoherence and concomitant fluidity 
observed in our numerical simulations implies thermalization of the system. We 
first investigate the behavior of the ensemble of initial conditions in the Poincare 
phase plane for the toy case of \miform background field and fluctuations. One 
sees that the initially localized trajectories spread around a close loop filling 
the phase-space as one would expect for the phase-space density of a micro- 
canonical ensemble. For the toy example considered, the ensemble average of 
the trace of the energy momentum tensor can be expressed at large times as 
the time average along a single trajectory in the Poincare phase plane. For the 
scale invariant (f)'^ theory, this average is zero with the consequence that the 
energy momentum tensor becomes traceless resulting in a single valued relation 
between the energy density and the pressure. 

Going beyond the toy example, for the general case of spatially non-uniform 
fluctuations, there is no easy way to visualize trajectories on the Poincare phase 
plane because the system is infinite dimensional. However, because the numer- 
ical problem is formulated on a lattice, one can look at a small sub-system on 
this lattice and study its event-by-event energy fluctuations. Starting from a 
Gaussian initial distribution, we see that the distribution converges to an expo- 
nential form. One can also study the moments of the energy distribution; these 
again demonstrate a rapid change from initial transient values to stationary 
values. While the behavior is close to those expected from a canonical thermal 
ensemble, it is premature from our present studies to make definitive conclu- 
sions. This will require a careful study of the effects of varying the coupling and 
volume effects and will be left to a future study. 

We note however that the formalism developed in our paper is well suited 
to the study of thermalization of quantum systems^. It has been argued pre- 
viously [38-41] that quantum systems will thermalize if they satisfy Berry's 
conjecture [42]. This conjecture states that the high lying quantum eigenstates 
of a system whose classical behavior is chaotic and ergodic have a wavefunc- 
tion that behaves as a linear superposition of plane waves whose coefficients 
are Gaussian random variables. When an inclusive measurement is performed 

^We thank Giorgio Torrieri for bringing to our attention Berry's conjecture and the ac- 
companying literature on eigenstate thermalization. 



4 



on such an cigcnstatc, one obtains results that agree with the predictions of 
the micro-canonical equilibrium ensemble, a property that has been dubbed 
"eigenstate thermalization" in [39] . If the state at t = is a coherent superposi- 
tion of such cigcnstatcs, the micro-canonical predictions become valid only after 
the states in the superposition have sufficiently decohered-thus for quantum 
systems where Berry's conjecture apply, thermalization appears to be a conse- 
quence of dccohcrcnce. The ensemble of quantum fluctuations included via the 
resummation we develop in the section 2.6 leads to fields that have precisely 
this behavior (see eqs. (51-52)). Our interest ultimately is in QCD, where the 
classical behavior of the system is believed to be chaotic [43 45] . Because much 
of our formalism can be extended to gauge theories, we anticipate that a first 
principles treatment of thermalization is feasible. 

This paper is organized as follows. In section 2, we introduce the model 
scalar theory and the CGC-like initial conditions for its temporal evolution. We 
then discuss the computation of T^" at leading and next-to-leading order. The 
problem of secular divergences is noted, and a stable resummation procedure is 
developed. A simplified toy model is considered in section 3, wherein only spa- 
tially uniform fluctuations are considered. The behavior of the resummed pres- 
sure and energy density and their relaxation to an equation of state is studied. 
These results are interpreted and understood as a consequence of the decoher- 
ence of the system which allows one to equate ensemble averages to a temporal 
average over individual classical trajectories. For the longitudinally expanding 
case, temporal evolution in the toy model displays the behavior of a fluid un- 
dergoing ideal hydrodynamic flow. The full quantum fleld theory is considered 
in section 4, where we compute ab initio the spectrum of fluctuations. The full 
theory displays the same essential features as the toy model studied in section 
3, albeit the interplay of linearly growing low lying momentum modes and the 
resonant modes leads to a more complex temporal evolution. In this section, 
we also investigate the dependence of the relaxation time on the strength of the 
coupling constant. Then, we study the energy distribution in a small subsystem, 
and its time evolution. We conclude with a brief outlook. Much of the details 
of the computation are given in appendices. In appendix A, we discuss the 
numerical solution of the scalar fleld model, including the lattice discretization, 
the computation of the quantum fluctuation spectrum and the sensitivity of the 
results to the ultraviolet cut-off. The stability analysis of linearized perturba- 
tions to the classical fleld is considered in Appendix B. The resonance band is 
identified and the Lyapunov exponents are computed explicitly. We also discuss 
the relationship between decoherence and linear instabilities. 

2 Temporal evolution of T^''' 

In this section, we will consider a scalar field toy model whose behavior mimics 
key features of the Glasma [24] description of the early behavior of the quark- 
gluon matter produced in high energy heavy ion collisions. In the CGC frame- 
work, strong color fields are present in the initial conditions for the evolution of 
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the Glasma. In this situation, the leading order contribution is given by classi- 
cal fields, with higher order corrections coming from the apparently sub-leading 
quantum fluctuations. We consider the stress-energy tensor in this scalar field 
model and discuss its temporal evolution at leading (LO) and next-to-leading 
orders (NLO). We show explicitly that there are contributions at NLO that can 
grow with time and become larger than the LO terms. We end this section by 
describing how these "secular" terms can be rcsummcd and the results expressed 
in terms of an average over a Gaussian ensemble of classical fields. 

2.1 Scalar model with CGC-like initial conditions 

Our CGC inspired scalar model has the Lagrangean 

/;^i(a^0)(5^<^)-^</.4 + j</>, (1) 

where J is an external source. In the CGC framework, the source J coupled to 
the gauge fields represents the color charge current carried by the two colliding 
heavy ions. The current is zero at positive proper time, corresponding to times 
after the collision has taken place. We emulate this feature of the CGC in 
a simpler coordinate system by taking the source J to be nonzero only for 
Cartesian time < 0, and parameterize it as^ 

j{x)^e{-x°)^. (2) 

At > 0, where J is zero, the fields evolve solely via their self-interactions, in 
an analogous fashion to the non-Abelian color fields produced in the collision of 
two hadrons or nuclei. 

In eq. (2), we incorporated two additional features of the CGC. The first 
feature corresponds to a strong external current J, which follows from the power 
of the inverse coupling when g <C 1; weak coupling is essential to motivate an 
expansion in powers of g^. The other feature of the CGC that is emulated is 
that the dimensionful parameter Q in eq. (2) plays a role analogous to that of 
the saturation scale [19,20], in the sense that non-linear interactions are sizeable 
for modes |fc| < Q. 

Note that a scalar field theory with a (j)'^ coupling in four space-time di- 
mensions is scale invariant at the classical level-the coupling constant g is di- 
mensionless in the theory. In our model, this scale invariance is broken by the 
coupling of the scalar field to the external source J containing the dimensionful 
scale Q. We may therefore anticipate that all physical quantities are simply 
expressed by the appropriate power of Q times a prefactor that depends on g. 

^In the numerical implementation of the model, the time dependent prefactor is constrained 
to vanish when —i- — oo to ensure a free theory in the remote past. 
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2.2 T^''' at leading order 



Because the source J contains a power of the inverse coupling, the power count- 
ing for Feynman diagrams indicates that the order of magnitude of a given graph 
depends only on its number of external lines and number of loops, but not on the 
number of sources J attached to the graph [46,47]. For the energy-momentum 
tensor of the theory, the various contributions can be organized in a series in 
powers of as 



rpfJtv _ ^ 
5^ 



Co + cig^ + C2g^ + ■ 



(3) 



In this expansion, the coefficients co,ci,C2, - - are themselves infinite series 
in the combination gj corresponding to an infinite set of Feynman diagrams. 
This combination is parametrically independent of g because J ~ g~^ . More 
precisely, Cq contains only tree diagrams, ci 1-loop diagrams, C2 2-loop diagrams, 
and so on. 

In our model, the leading order (tree level) contribution to the energy- 
momentum tensor can be expressed solely in terms of a classical solution ^p 
of the field equation of motion [46]. Namely, one has 



where 



(4) 



lim <^(a; , cc) = . 



(5) 



Clearly, due to the non-linear term in the equation of motion, the solution 

(and hence the coefficient cq) depends on gJ to all orders, as stated previously. 
This LO energy momentum tensor is conserved^. 



3 T^" 



. 



(6) 



If the source J is taken to be spatially homogeneous, then the energy- 
momentum tensor evaluated at leading order has the simple form 

) 





V 

with the leading order energy density and pressure given by 

1 



(7) 



2^ 4!^ 

1-2 5^ 4 

= 2^ -4!^ 



(8) 



^Strictly speaJcing, this is true only at > 0. At negative times, some energy is injected 
into the system by the external source J . 
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Figure 1: Components of Tj^^ for a spatially uniform external source. To per- 
form this calculation, we took in eq. (5) a source J = g~^Q^9{—x^)e'''^^ (with 
5=1,6 = 0.1 and Q ~ 2.5), that vanishes adiabatically in the remote past. 



One can easily check that the energy density e^^^ is constant in time at x'' > 
(after the external source J has been switched off), while the pressure p^o 
a periodic function of time at > 0, as illustrated in the figure 1. From the 
numerical computation, it is clear that at this order of the calculation of e^^^ 
and Pj^Q , one does not have a well defined (single valued) relationship e^^^ = 
/(Plo)- other words, there is no equation of state at leading order in g^. 
This might appear problematic at the outset because one might expect that the 
scale invariance of the theory would require the energy momentum tensor to be 
traceless. As discussed further in section 3.4.2, this is not so for the case of a 
scalar theory. 

2.3 T^" at next to leading order 

At next-to-leading order, the energy momentum tensor can be written as 
(fik 



J (27r)32fc —ydaa-kd"a+k ~ V"{ifi)a-ka+k 



(9) 



where for brevity we use the notation V{(p) = g'^ip^ /Al with each prime denoting 
a derivative with respect to f. In this formula, /? and a±k are small field 
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perturbations, that are defined by the following equations: 



□ + V"{v) 



a±k = 



lim a±k{x) = e 



lim P{x) = . 



(10) 



Because the classical field (f is spatially homogeneous in the toy model consid- 
ered here, the equation of motion for a±k simplifies to 



a±k + (fe' + V"iip))a±k = 



(11) 



and the field fiuctuation j3 depends only on time. 

After some algebra, it is easy to check that the energy-momentum tensor is 
also conserved at NLO^ for > 0, 



(12) 



The 00 component of T^" in eq. (9) gives us the energy density at NLO, 



^0 + /3y'(^) + - 



2 J (27r)32fc 



a-kCb+k + (fe^ + V" {ip))a-ka+k ■ (13) 



Given eqs. (10), it is straightforward to verify that this correction is also constant 
in time, e^Lo = 0> agreement with eq. (12). The 11 component of eq. (9) -the 
NLO pressure in the x direction- reads 



2 J (27r)32fc 



-ka+k - (fc^ - 2fc^ + V"{ip))a-ka+k 



(14) 

Note that although the integrand is not rotationally invariant, the result of the 
k integration is symmetric and the NLO pressures are the same in all directions. 

We evaluated numerically ej^^o '^^^'^ ^'nlo ^ coupling constant 5 = 1, by 
first solving eqs. (10) for /9 and for the a^'s (for a discretized set of fe's). The 
results of this calculation are shown in the figure 2. From this evaluation, we 
see that the energy density at NLO is constant at x° > 0, as we expected^. 
We also notice that for g = 1, the NLO correction to the energy density is 
very small, of the order of 1.4% of the LO result^. Thus, we conclude from 



^This result should be self-evident because the conservation equation djj.T'^'^ = is linear 
in the components of T*^". Therefore, it does not mix the different orders, requiring the 
conservation equation to be satisfied for each order in p-^. 

®This time independence can be seen as a test of the accuracy of the numerical calculation, 
because it results from a cancellation between several terms that grow with time. 

^Indeed, since there is a prefa<;tor 1/4! in our definition of the interaction potential, g = I 
corresponds to fairly weak interactions. 
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Figure 2: Components of T^'^o ^'^^ ^ spatially uniform external source. This 
calculation was performed for = 1. 



this that for such a value of the coupling, we have a well behaved perturbative 
expansion for e. The NLO pressure however behaves quite differently. Not 
only it is varying in time (hence no equation of state at NLO), but it also 
has oscillations whose amplitude grows exponentially at large x^. Therefore, 
the NLO correction to the pressure eventually becomes larger than the LO 
contribution, and the perturbative expansion for the pressure in powers of 
breaks down^. Also noteworthy is the fact that at — 0, p^^^ is still a small 
correction to p^o ! only becomes large at later times. 

2.4 Interpretation of the NLO result 

The secular divergence of the pressure at NLO can be understood as a conse- 
quence of the unstable behavior of a±;,(x) for some values of k. The stability 
analysis of small quantum fluctuations in 0^ field theory is performed in ap- 
pendix B. From this study, one obtains the following results: 

i. There is a range in |fc| where the aj-^'s diverge exponentially in time, due 
to the phenomenon of parametric resonance. 

ii. The zero mode fe = fluctuation, ao, diverges linearly in time, a phe- 
nomenon closely related to the fact that the oscillation frequency in a 
non-harmonic potential depends on the amplitude of the oscillations. 

In addition, one observes numerically that fluctuation modes in the vicinity of 
k — 0, albeit not mathematically unstable, can attain quite large values. (They 
appear to grow linearly for some time before decreasing in value.) 

similar behavior was observed in a different context in [48] . 
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Because of the existence of modes that grow in time, integrals such as 

/(fik 
^2^)32fc a-k{x)a+k{x) , (15) 

that appear in the components of T^^^ (see eqs. (13) and (14)) or in the right 
hand side of the equation (eq. (10)) for /3, are divergent when — >■ +oo as 
illustrated in the figure 3. In this plot, one can check that the envelope of the 
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exp(lt) 1(11 

Figure 3: Numerical evaluation of the integral defined in eq. (15). The line 
denotes an exponential fit to the envelope. 

oscillations grows exponentially, with a growth rate A « 2 * //max where /imax 
is the maximal Lyapunov exponent in the resonance band. If the integral in 
eq. (15) is evaluated with an upper cutoff that excludes the resonance band 
from the integration domain, then I{x'^) grows only linearly, because now its 
behavior is dominated by the soft fluctuation modes whose growth is linear. 

Even though secular divergences in integrals such as eq. (15) are present in 
eq. (13), they cancel in the calculation of e„Lo because the energy density in 
our toy model is protected by the conservation of the energy momentum tensor. 
However, they do not cancel in p^lo which explains the divergent behavior 
displayed in fig. 2. 

2.5 Alternate form of T^" at NLO 

The secular divergence of the pressure at NLO suggests that the weak coupling 
series for the pressure may be better behaved if one develops a resummation 
scheme that captures the physics of the secular terms by identifying their con- 
tribution and summing them to all orders in perturbation theory. Before we do 
this, we shall discuss a general formulation of the energy-momentum tensor at 
NLO which will help formulate the problem of resumming secular terms. 
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In previous works [46,47], wc showed that the problem of computing NLO 
corrections for inclusive quantities-such as components of the energy momentum 
tensor in field theories with strong sources could be formulated as an initial 
value problem. Specifically, for the cncrgy-momcntiim tensor, wc can write the 
NLO contribution at an arbitrary space-time point as the action of a functional 
operator acting on the LO contribution, 

(16) 

The operator that appears in eq. (16) is the generator of shifts of the initial 
conditions ifiojdofo (at x° = 0) of the classical field, 

S S 
a-T^ = a{0,u)——- + d{0,u)— ^. (17) 

The factor in the functional formulation of eq. (16) should therefore be 
considered as a functional of the value of (p, 93 at = 0. The full content of the 
temporal NLO evolution of T^'' is contained in eq. (16). One can check that 
this expression is exactly equivalent to eq. (9) [25]. 

The expression in eq. (16) has been obtained by splitting the time evolution 
at x'^ = such that the x'^ < part of the time evolution is described by the 
operator in the square brackets, and the evolution at x''^ > is hidden in the 
functional dependence of T^^^ with respect to the value of the classical field Lp 
at x'^ = 0. The choice of a;*' = for this split in the time evolution is arbitrary 
and equivalent formulas can be obtained with other choices^. Here, our choice 
is motivated by the fact that x'^ is the time at which the external source J turns 
off. In view of the rcsummation we will use later, it is important to note that 
the quantum field fluctuations /3 and a±k are still small relative to the classical 
fic^ld at the splitting time used in the formula. That this is true in our case is 
transparent from the figure 2. 

2.6 Resummation of the NLO corrections 

As seen previously, the fixed order NLO calculation is not meaningful after a 
certain time, because it gives a pressure that is larger than the LO contribution. 
The NLO contribution (and likely any higher fixed loop order contribution) has 
secular divergences because it involves the linearized equation of motion for 
perturbations to the classical field (p. In other words, if ^ = + o, the NLO 
calculation approximates the dynamics of by 

nip + v'{ip) = J 

a + V"{ip)\a = 0, (18) 



*The splitting of the time evolution in two halves need not be done at a constant and 
any locally space-like hypersurface will suffice. 
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on the grounds that the nonlinear terms in a arc formaUy of higher order in g^. 
Obviously, if the dynamics of tp was treated exactly, by solving instead^ 



+ V'{il)) = J , (19) 

we would not have any divergence because the tp^ potential would prevent run- 
away growth of tp- However, in order to achieve this substitution, we must 
include in our calculation some contributions that are of higher order in g^. 
Thus, wo seek a rosummation that restores the lost nonlincarity in the field 
fluctuations, while keeping in full the LO and NLO contributions that we have 
already calculated. 

As we will argiie in this section, a simple rcsummation that leads to an 
energy-momentum tensor which is finite at all times consists in starting from 
eq. (16) and in exponentiating the operator inside the square brackets, 

TIZum{x) = exp[ jd'up-Tu+\ j d^ud^v j ^^^!^[a+fe-T„][a_fe-T„]]T^';^(a;) , 

(20) 

If we Taylor expand the exponential, we recover the full expressions for the LO 
and NLO contributions, plus an infinite series of other terms that are of higher 
order in g'^ , 

TZu..{x) = % [ CO + ciff2 + + •••]. (21) 
fully partly 

From the form of eq. (20), it is not evident that the exponentiation leads to a 
better behaved result; on the surface it appears that we are including an infinite 
series of terms that are increasingly pathological at large times. To see that the 
result is now stable when ~> +00, let us consider some generic function of 
the classical field at the point a;, F[ip{x)]. The field (p{x) is itself a functional of 
the values^° ipo of the field and (fo of its first time derivative at = 0. Thus, 
the quantity F[(fi{x)] is implicitly a function of ipQ,ipQ, 

F[<p{x)]=F[<po,V^o]. (22) 

Note now that the exponential of 13 ■ T„ is a translation operator when it acts 
on a functional F[(pQ{u),(fio{u)], 



exp 



j d^uP- T„] F{^o, ^0] = F[^o + P,^o+P]- (23) 



The first term in the exponential in eq. (20) therefore merely shifts the initial 
conditions ipo, (po at x*^ = of the classical field ip (by amounts /3). Similarly, 



^Though this expression looks identical to the first equation of eq. (18), the initial condi- 
tions for this equation are different, leading to a difl'erent solution. 

^''Although we do not write that explicitly in order to simplify the notation, ipo and ipo may 
depend on the position x. 
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the second term, that involves the exponential of an operator that has two T's, 
can be rewritten as a sum over fluctuations of the initial classical field"'^"'^ 



exp 



If f (Pk 1 

- j cPud'v j j^-^[a+k ■ T„][a_fe • T„]J F[ipo,^o] = 

= J [DaDa] Z[a, a] F[<fo + a,<po + a], (24) 

where the distribution Z[a,a] is Gaussian in a (a;) and a{x), with 2-point cor- 
relations given by 

f (fik 

{a{x)a{v)) = J j^--^a+k{0,x)a-k{0,y) , 

/(Pk 
^^-^a+k{0,x)d-k{0,y) . (25) 

Therefore, the energy-momentum tensor resulting from the resummation of 
eq. (20) can be written as 

nZun. = J [Da{x)Da{x)] Z[a, a] T^^[v^ + /3 + a] , (26) 

where T^^ [i^o + /3 + ct] denotes the LO energy-momentum tensor evaluated with 
a classical Held whose initial condition at x° = is + /3 + a (and likewise for 
the first time derivative). 

From eq. (26), one can now see why the proposed resummation cures the 
pathologies of the NLO contribution. While the fixed-order NLO result involved 
linearized perturbations to the classical fields (that are gcnerically divergent 
when x° oo), in the resummed expression these perturbations appear only as 
a shift of the initial condition for the full non-linear equation of motion. After 
this resummation, the evolution of the perturbations at > is no longer 
linear-since the (j)'^ potential is bounded from below the evolution is stable. 

In addition to manifestly demonstrating the stable evolution demanded by 
the undcaiying theory, eq. (26) is a most useful expression for a practical im- 
plementation of our resummation. It is important to note however that the 
integral over fc in the 2-point correlations (eq. (25)) that define the Gaussian 
distribution of a and a should be cut-off at a value A ~ gipo ^ Q in order 
to avoid ultraviolet singularities. With such a cutoff, one can show that the 
sensitivity to the value of the cutoff is of higher order in g'^, while at the same 
time being large enough to include in the resummation all the relevant unstable 
modes (the modes with Q ^ \k\ are all stable). 

^^An elementary form of the identity, 

-,2 e-^^/^f 
e2 9./(x)=/ dz fix + z), 

can be proven by doing a Taylor expansion of the exponential in the left hand side and of 
f{x H- z) in the right hand side. In this simple example, one sees that an operator which is 
Gaussian in derivatives is a smearing operator that amounts to convoluting the target function 
with a Gaussian. Another way of proving the formula is to apply a Fourier transform to both 
sides of the equation. 
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3 T^esum from spatially uniform fluctuations 



Before we proceed to a full 3+1-diniensional numerical evaluation of eq. (26) 
with an ab initio computation of cq. (25), wc shall first consider, as a warm- 
up exercise, a computation including only spatially homogeneous fluctuations. 
Albeit not realistic, this much simpler calculation will be very instructive in 
understanding the effects of these fluctuations on the behavior of the energy- 
momentum tensor. 



3.1 Setup of the problem 

For spatially homogeneous fluctuations, the main simpliflcation is that func- 
tional integrations over the fields a and a in eq. (26) become ordinary integrals 
over a pair of real numbers, with the Gaussian weight 



Z{a, a) = exp 



2(Ji 2(j2 



(27) 



The two parameters C7i,2 can be used in this toy calculation to control the 
magnitude of the fluctuations. In the limit (Ji.2 0, we recover the leading 
order result which of course receives no contribution from the fluctuations. 

The second important simplification in this toy calculation is that since both 
the underlying classical field and the fiuctuations are spatially homogeneous, 
the field equation of motion is an ordinary differential equation^^. One should 
note that the characteristic oscillation frequency is directly proportional to the 
amplitude of (p{t). This property of the solution will be key in interpreting the 
results that follow. 

We now turn to the computation of the resummed pressure and energy 
density in this toy model. Prom eqs. (8) and (26), the expressions for the 
energy density and the pressure read 

lv>'-Vi^)) , (28) 



^^In this case, the field equations can even be solved analytically. This can be seen very 
simply: from energy conservation, i (p"^ + V{<p) = Eq = y(^max) (with (/)max the amplitude of 
the oscillations of ip{t)), on gets 



- const + 



V2J0 



V2J0 sfvT^ max 

) - 



For a potential, the integral in the right hand side is an elliptic integral, and one can 

express as 

ipit) = fimax CIli/2 (9</'maxt/\/24 + COUSt) , 

where cni/2 is the Jacobi elliptic function of the first kind with the elliptic modulus k = 1/2. 
This expression is periodic with a period T = 2\/24i<'(l/2)/g(^inax, where i<r(l/2) « 1.85 is 
the complete elliptic integral of the first kind. 
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where ip is the solution of the classical equation of motion whose value at = 
is (/3o + a and whose time derivative at = is 930 + The brackets (•••)„ 
denote an averaging over all possible values of a, a with the distribution of 
eq. (27). 

3.2 Energy momentum tensor 

In fig. 4, we display the result of the toy model calculation in the limit where 
we do not have fluctuations (0-1^2 0). As anticipated, the result is equivalent 
to the one displayed in fig. 1 for the leading order calculation. In this figure, 




20 40 60 80 100 

time 



t" T»»/3 

Figure 4: Components of Tlq, where when no quantum fluctuations are in- 
cluded. 

for reasons that will become obvious shortly, we have represented the energy 
density divided by three. In fig. 5, we show the results of the same calculation 
performed with non-zero widths ai^2 for the Gaussian distribution of fiuctua- 
tions. We observe a striking difference of the resummed result compared to the 
previous (LO) figure-the oscillations of the pressure are damped and the value 
of the pressure relaxes to e/3. Subsequently, one has a single- valued relation- 
ship between the pressure and the energy density, namely, an equation of state 
- specifically, the equation of state e = 3p of a scale invariant system in 1 -I- 3 
dimensions. 

3.3 Phase-space density 

It is also instructive to look at the phase-space density pt {(fi, (fi) of the points 
(if, ip) as the system evolves in time. This is shown in fig. 6. At t = 0, we start 
with a Gaussian distribution of the initial conditions, with a small dispersion 
around the average values {(p = 10 and <y9 = in our example). 



16 




Figure 5: Components of TJ-osum obtained with a Gaussian ensemble of spatially 
uniform quantum fluctuations. 



t = 200 t= 50 t- 10 Q t= 5 t- 



Figure 6: Phase-space distribution of the ensemble of classical fields at various 
stages of the time evolution. 
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Each initial condition then evolves independently according to the classical 
equation of motion, and the corresponding trajectory in the (tp, ip) plane is a 
closed loop^^ due to the periodicity of classical solutions. One observes that the 
initially Gaussian-shaped cloud of points starts spreading around a closed loop, 
to eventually fill it entirely when a;° — +00. When this asymptotic regime is 
reached, the density pt{'f, 0) depends only on the energy-roughly speaking, the 
radial coordinate in the plot of figure 6 and no longer on the angular coordinate. 

A more formal way of phrasing the same result is to first note that the time 
evolution of the phase-space density pt obeys the Liouville equation, 

^ + {p„iJ} = 0, (29) 

where {•, •} is the classical Poisson bracket. Therefore, if a stationary distribu- 
tion is reached at late times, it can only depend on (p and (p via H{ip^ ip). The 
asymptotic behavior of the phase-space density in our toy model is reminiscent 
of a micro-canonical equilibrium state, in which the phase-space density is uni- 
form on a constant energy manifold^^. In other words, all micro-states that 
have the same energy are equally likely. 

3.4 Interpretation of the results 

We shall now discuss the physical interpretation of our results, first discussing 
the decoherence of the temporal evolution of the fields and their time derivatives, 
and subsequently, the impact of decoherence on the relaxation of the pressure 
towards that of a scale invariant system. 

3.4.1 Decoherence time 

Of the previous numerical observations, the easiest to understand is the spread- 
ing of the phase-space density around a closed orbit. Because the oscillations 
are non-harmonic, the various points in the plot of figure 6 rotate at different 
speeds^^; in a </)^ potential, the outer points rotate faster than the inner ones. 
Therefore, as time increases, the cloud of points spreads more and more due to 
this effect. 

One can estimate the time necessary for the cloud of points to spread over a 
complete orbit. This happens when the angular spread of the points reaches the 
value 27r. For one field configuration, this angular variable is, up to a phase that 
depends on the initial condition, 6 — tut, and the angular velocity to depends 

-'^These loops arc constant energy curves ^ip'^ + V{ifi) = H. 

^"^It should be noted here that a spatially homogeneous field is very special regarding this 
issue; indeed, any non-linear system with a single degree of freedom is ergodic. This is not 
necessarily the case if there are more than one degrees of freedom, as is the case in a full 
fledged field theory. 

-'^^Thc assumption of a scale invariant theory simplifies some expressions here, but is not 
crucial to the argument. The only requirement for this phenomenon is that the frequency 
of the oscillations depends on their amplitude; thus any non-harmonic potential will lead to 
similar results. 
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only on the energy of that particular field configuration. (In our case, this 
phase is small for a narrow Gaussian distribution.) If we consider two field 
configurations, their angular variable difference increases linearly in time, 
A0 = Auj t, where Alj is the difference between their angular velocities. In the 
case of a 5^(/)^/4! potential, one can prove that (see the footnote 12) 

^=:r7^ Jf"7, ~ 0-346 5<^max , (30) 



VT- 



where 0max is amplitude of the oscillations of the ip field. Thus, the angular 
shift between the two field configurations is also w 0.346 gA<pjnaxt, and this 
shift reaches 2tt in a time 

18.2 , , 

t^——. (31) 

.gA^max 

After this time, the two fields have become completely incoherent. We see that 
this time is inversely proportional to the coupling constant g, and to the differ- 
ence of the field amplitudes. Thus a narrow initial Gaussian distribution will 
need a longer time to spread around the orbit than a broader initial distribution. 

3.4.2 Equation of state from quantum averaging 

Once we know that the phase-space density spreads uniformly on constant en- 
ergy curves, it is easy to understand why the pressure relaxes towards e/3 when 
we let the initial conditions for the classical field fluctuate. The trace of the 
energy-momentum tensor (assuming 4 dimensions of space-time) is 

T% = <^(^n<^ + 4^^ -5c.((p5"v5) . (32) 

A scale invariant theory in four dimensions is a theory in which the interaction 
potential obeys V'{(p) = AV{ip)/ip. This is the case of a (/)^ interaction. There- 
fore, the first term in the right hand side of the previous equation vanishes 
thanks to the equation of motion of the classical field ip. This result shows that 
the energy-momentum tensor of a single configuration of classical field is not 
zero in our model, but is a total derivative-'^®. In our simplified toy model where 
the fields are spatially homogeneous, the previous relation simplifies to 

T% = -«. (33) 

When averaged over one period, the trace of the energy-momentum of one clas- 
sical field configuration vanishes because the classical field is a periodic function 
of time, 

T^.^^l dTT%(^(T),^(r))=0, (34) 



^^There is cui alternative "improved" definition of the energy-momentum tensor that is 
explicitly traceless [49]. However, while the energy density has a single valued relation to 
the pressure, this pressure is not the canonical pressure. As we shall discuss later in section 
3.5.3, both definitions give a deviation from ideal hydrodynamic flow, which is cured by the 
quantum averaging described here. 
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where the result is independent of t. When we calculate the energy-momentum 
tensor averaged over fluctuations of the initial conditions, we are in fact per- 
forming an ensemble average weighted by the phase-space density pti^p, <p), 

(^%)a,d = pM^) T%{^,^) , (35) 

and the time dependence of the left hand side comes from that of the density pt- 
It is convenient to trade the integration variables (p, if for energy/angle variables 

E,e, 

= J dEde ptiE, e) T%{E, 9) , (36) 

where pt is the phase-space density in the new system of coordinates^^. Our 
first result shows that pt{E, 9) — > pt{E), namely, becomes independent of 9 at 
late times, which enables us to write the previous equation as 

(^%)a,d , / dE Pt{E) J d9 T^,{E, 9) . (37) 

The crucial point here is that the integral over 9 is simply the integral over one 
orbit for a single classical field configuration (eq. (34)), 

/o pt-\-T' 
dO T%{E, 0) = yI '^'^ T\{^{t), ^{t)) = . (38) 

Thus, we have proven that 

e-3p=(T%)„_. ^«^0, (39) 

in agreement with what we have observed numerically. Moreover, from the 
derivation of this result, it is clear that the time necessary to reach this limit 
is the same as the time (in eq. (31)) necessary for the phase-space density to 
become independent of the angular variable 9. 

3.5 Effect of the longitudinal expansion 

We have thus far considered a system of strong fields enclosed in a box of 
fixed volume. There is therefore no concept of hydrodynamical flow in such a 
system. To fully understand the implications of decoherence and relaxation of 
the pressure we have discussed previously for hydrodynamical flow, we will now 
simply generalize the toy problem of spatially uniform fields and fluctuations to 
a system undergoing a boost invariant one dimensional expansion. 

^"^ pt is equal to the original pt times the Jacobian of the change of variables. 



20 



3.5.1 Relaxation of the pressure 

The geometry of the one dimensional expansion (chosen to be the z direetion) is 
appropriate to describe the coUision of two projectiles (nuclei) at ultrarelativistic 
energies. The natural coordinates are the proper time r and rapidity r] defined 
by 
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, . iln(l±f). (40) 

In the spatial plane orthogonal to the z axis, the coordinates are denoted by 
x^. In this system of coordinates, the classical equation of motion for a field y 
that depends only on proper time is 

^+-<^+^<^^ = 0, (41) 

T 

where the dot now denotes a derivative with respect to r. The analog of the 

the toy problem we discussed previously in the first part of this section is to 
let the initial conditions of the field have Gaussian fluctuations a, d that are 
independent of r] and x^. 

The components of the energy-momentum tensor in this system of coordi- 
nates, averaged over the fiuctuations of the initial conditions are 

\ 2 I a, a 

As in the fixed volume case, we shall use the distribution of eq. (27) for a, a. 
The result of this computation is shown in the figure 7. The dots represent 
the energy density divided by 3, and one observes that its time dependence is 
well described by a r~^/^ decay characteristic of boost invariant flow in ideal 
relativistic hydrodynamics. If we do not include fluctuations of the initial con- 
ditions, we observe that the pressure oscillates between positive and negative 
values, with a decreasing envelope. Conversely, if we average over an ensem- 
ble of initial conditions, we see the oscillations of the pressure dampen quickly, 
ensuring that the pressure approaches one third of the energy density. 

These results are in sharp contrast to what one obtains for a (fy^ potential. 
The results are shown in flg. 8. In this case, the fluctuations do not make the 
pressure converge to T°°/3, and the latter decreases as instead of r""^/^. 

3.5.2 Interpretation of the results for expanding fields 

From the equation of motion in eq. (41), we obtain 
de 



e + p 



1 .2 

--(p 

T 



(43) 
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{|)* potential, longitudinal expansion 
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|T^| w/o fluctuations 

|T^| with fluctuations 
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Figure 7: Numerical evaluation of T^^ for fields undergoing a boost invariant 
1-dimensional expansion in a (j)^ theory with a Gaussian ensemble of spatially 
uniform initial fluctuations. 



(|)^ potential, longitudinal expansion 
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Figure 8: Numerical evaluation of T^'^ for fields undergoing a boost invariant 
l-dimensional expansion in a (jP theory, with a Gaussian ensemble of uniform 
initial fluctuations. 
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This equation, which is vahd for individual classical field configurations at every 
time, is identical to Euler's equation for boost-invariant ideal hydrodynamics. 
The difference with hydrodynamics lies in the fact that hydrodynamics assumes 
the existence of an equation of state p = f{e) to ensure a closed form expression 
in eq (43). In classical field dynamics, one is not free to impose a relationship 
between e and p since they are both completely determined from the field ip and 
its derivative ip. For instance, as seen in fig. 7. for a single classical solution, we 
do not have a one-to-one correspondence between e and p; e has a monotonous 
behavior while p oscillates. What is remarkable is that the ensemble average 
over the initial conditions leads in a short time to a one-to-one correspondence 
e = 3p, which is precisely the equation of state one would use in boost invariant 
hydrodynamics of a perfect fluid. 

The mechanism whereby this relationship is reached is the same as in the 
non-expanding case. As previously, one can prove that for a single phase-space 
trajectory, the time averages of e and p obey a relation identical to the expected 
equation of state, because the trace of the energy momentum tensor is a total 
derivative. Then, by using the fact that different initial conditions lead to 
different oscillation frequencies, one gets the phase decoherence that enables us 
to transform the ensemble average over the initial conditions into a time average 
along one classical field trajectory. This decoherence is the missing ingredient 
in the harmonic case-as we noted previously, it arises for the theory because 
the angular velocity of the phase space trajectory of an individual conflguration 
depends on the amplitude of the configuration. 

From this result, it is very easy to obtain the t~^/^ behavior of the energy 
density. The ensemble average of eq. (43) at late times is 

which leads immediately to the observed behavior. Since both e and p decrease 
like T^^/^ even for a single configuration (if one considers the envelope of the 
oscillations of p), this means that at late times we have 

if ~ , if ~ r"^/^ . (45) 

This behavior is seen from the simple ansatz (/?(r) ~ cos(/(r))T~-^/'^, which, 
while inaccurate in detail, qualitatively captures the right physics. For a 
potential, the frequency f (x (p ^ t--i/3 for a (j)^ potential; using this relation in 
our ansatz gives the stated result. The fact that cp decreases faster than f is due 
to the slowing down of the oscillations with time, as their amplitude decreases. 



3.5.3 Callan- Coleman- Jackiw energy-momentum tensor 

An alternate definition of the energy-momentum tensor was proposed by Callan, 
Coleman and Jackiw (CCJ) [49]. Their expression is explicitly traceless. They 
argued further that their form of the stress energy tensor improved properties 
relative to the canonical one with regard to renormalization. Let us briefly 
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summarize the differences between the usual definition of T^'^ and CCJ's. With 
the canonical definition that we have used thus far, one has, 



dr 

de e + p 



(46) 



dr T 

With this definition of T'"', one obtains the equation for Bjorken hydrodynamics 
automatically for each configuration of the classical field. However, one gets 
e = 3p only through decoherence, by averaging over an ensemble of initial 
conditions. 

In comparison, with CCJ's definition of the energy-momentum tensor, one 

has 

P = \^^-V{^)-'^Uip^ 

e - 3p = 

de ^ e+p 1 djifiip) 

dr T 3t dr ' ^ ' 

With this form of the energy momentum tensor, the equation of state is satisfied 
for each classical field configuration, but not Bjorken's hydrodynamic equation. 
It is only after an average over an ensemble of initial conditions that the last 

term d{(fi\p) /dr) in the last equation vanishes by decoherence. 
Since one requires simultaneously 

^ + 1±P = 
dr T 

e-3p = , (48) 

for ideal hydrodynamical flow, one sees that there is no discrepancy between 
the two descriptions despite the apparent differences. For our choice of the en- 
ergy momentum tensor, the reason why we didn't have an ideal hydrodynamic 
behavior at the beginning of the evolution was because of the lack of an equa- 
tion of state. In the case of CCJ's energy momentum tensor, it is because of a 
violation of the canonical Euler equation. The net effect of the quantum aver- 
aging in each case is to get rid of one or the other violation thereby ensuring 
ideal hydrodynamical behavior. In the following section, we will discuss only 
the canonical energy momentum tensor, for which the focus is on obtaining the 
equation of state as the necessary condition for hydrodynamical flow. 
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4 Results from the full fluctuation spectrum 



In the previous section, we showed that averaging over an ensemble of initial 
conditions for classical fields can lead the pressure to relax towards one third of 
the energy density. However, this study was oversimplified since we used only 
fluctuations that are uniform in space, and their Gaussian distribution was set 
by hand. However, quantum field theory predicts what the spectrum of these 
fluctuations is: one should average the LO energy- momentum tensor^^, 

TZum = {Ti:;:[^o + a])^ . , (49) 

over space-dependent random Gaussian fields a and d that have the following 
variance: 

which leaves no freedom to handpick what fluctuations we use. Prom these for- 
mulas, we can numerically compute ab initio the behavior of the pressure. The 
only tunable quantities in the calculation are then the scale Q (or more gener- 
ally the source J) that controls the amount of energy injected into the system 
at t < 0, and the coupling constant g. Note that the above spectrum of fluctua- 
tions for the initial condition for the field (p is equivalent to parameterizing the 
initial field as^^ 

ip{0,x) = (po{0,x) + / „„, Cka+k{0,x) + cla-k{0,x) , (51) 



(27r)32fc 

where the Cfe are random Gaussian numbers with the following variance 

(cfect)=0, {ckC*i) = {2Trf\k\S{k-l). (52) 

Details of the numerical lattice computation are relegated to the appendix A; 
in this section we focus on the results from numerical simulations with this ab 
initio spectrum of fluctuations. 



4.1 Numerical results 

Unless stated otherwise, the numerical results in this section are obtained on a 
12^ lattice^". The functional integration in eq. (26) is approximated by a Monte- 
Carlo average over 1000 configurations of the initial conditions, distributed ac- 
cording to eqs. (51) and (52). 

18 Since, at x" = 0, 13 (see eq. (26)) is a small shift that does not fluctuate, we have absorbed 
it into a redefinition of the classical field (po- 

i^If we recall that the a±k's are plane waves modified by the presence of the background 
field (fio, we observe that the fluctuating part of the initial fleld is very similar to the form of 
the wavefunction of high lying eigenstates for quantum systems that obey Berry's conjecture. 

^"In some instances, we have also performed simulations on a 20^ lattice and found only 
very small differences as long as the physical scales are below the lattice cutoff. 



25 



In fig. 9, we show the result of the computation of the pressure averaged 
over the Gaussian ensemble of initial conditions, for a value of the coupling^^ 
g ~ 0.5. We also show the energy density divided by three on the same plot. 
All the quantities in this plot are expressed in lattice units, which means that 




50 100 150 200 250 



Figure 9: Time evolution of the pressure averaged over the initial fluctuations. 
All the resonant modes are included in the simulation. The coupling constant 
is g = 0.5. 

the horizontal axis is t/a (where a is the lattice spacing) and the vertical axis 
should be understood as ea*/3 or pa'^. The lattice cutoff in this simulation is 
chosen to be just above the upper limit of the parametric resonance window 
(fc/mo = 3~^/^ where TOq = g'^ip^/2); therefore, all the resonant modes take 
part in the dynamics of the system. 

We observe that the ensemble averaged pressure relaxes towards e/3. This 
plot, obtained with the spectrum of fluctuations predicted by quantum field 
theory, is one of the central results of this paper. One can qualitatively identify 
two stages in this relaxation: (1) in the range < t < 50, the amplitude of the 
pressure oscillations decreases very quickly to a moderate value and, (2) from 
time 50 onwards, one has a slower approach of the pressure to e/3 that gets 
slowly rid of the residual oscillations. We will observe again this two-stage time 
evolution when we look at the fluctuations of the energy density. 

4.2 Influence of the resonant modes 

In section 3, we observed that the pressure relaxes to e/3 even if only the mode 
A; = is included in the simulation. This was understood as an consequence of 
the phase decoherence that exists in a non-harmonic potential between classi- 
cal solutions that have slightly different amplitudes. When we include all the 

^^Since the prefactor in the interaction potential is 3^/4!, a value g = 0.5 corresponds to a 
very wreak coupling strength. 
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Figure 10: Time evolution of the pressure averaged over the initial fluctuations. 
The lattice cutoff is located below the resonance band in order to exclude them 
from the simulation. The coupling constant is g — 0.5. 



fc-modes of the fluctuations, the situation becomes more complicated. In partic- 
ular, the stability analysis of these fluctuations (see the appendix B) indicates 
that in addition to a linear instability of the soft modes due to the above men- 
tioned decoherence phenomenon, there are also exponentially unstable modes 
in a narrow band of values fc. 

In order to assess the role played in the time evolution by the modes of 
the resonance band, we performed a second simulation with the same physical 
parameters, but now with the lattice cutoff placed just below the lower end of 
the resonance band. This makes certain that none of the modes that exist on 
this lattice has an exponential instability. Since the resonance band is quite 
narrow, this is a small change of the cutoff in physical units because the cutoff 
in the earlier simulation was just above the upper end of the resonance band. 
However, one can see in the figure 10 that excluding the resonant modes leads 
to significant changes. 

The final outcome, the relaxation of the pressure towards e/3, is not changed, 
but the details of the time evolution of the pressure are modified. Firstly, one 
observes a rather long delay during which the oscillations of the pressure remain 
almost constant in amplitude. Then, at a time of order 75 in lattice units, 
these oscillations are damped very quickly to very small wiggles around e/3. 
Except for a brief relapse, the oscillations remain very small after this time. In 
particular, the two-stage evolution that we observed with the full spectrum is 
now replaced by the following two stages: (1) nothing happens and, (2) very 
rapid relaxation that leaves almost no residual oscillations. 

Therefore, it appears that the resonant modes, even if their presence or 
absence in the resummation does not change the final outcome, do alter signif- 
icantly the detailed time evolution of the pressure. At this point, the precise 
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role of the resonant modes is somewhat unclear. It appears that the dynamics 
of the complete system is much richer than what one can learn by studying the 
linearized evolution of a single mode as done in the stability analysis of the ap- 
pendix B. This analysis does not capture the non- linear couplings between the 
various modes (once the instabilities have made them large) which gives them a 
big role in the late stage evolution of the system. This certainly deserves further 
study. 

4.3 Dependence on the coupling constant 

The simulation that led to the result of fig. 9 was performed with a value g = 0.5 
for the coupling constant - a very small value for our scalar field theory since 
there is also a factor 1/4! in the interaction potential. We have studied the time 




50 100 150 200 250 

time 

9=0.5 g=1 g=2 g=4 g-8 e/3 

Figure 11: Time evolution of the pressure averaged over the initial fluctuations 
for various values of the coupling constant: g = 0.5,1,2,4,8. All the resonant 
modes are included in the simulation. See footnote 22. 

evolution of the pressure for larger values of the coupling constant: g — 1,2,4,8, 
and the results are shown in the figure 11. Note that this computation is done 
at fixed energy density. Indeed, because Q is the only dimensionful parameter of 
our model and there is a factor 1 /g in the source J, the energy density behaves 
at leading order as e cx Q'^/g^. Thus, if we increase g at constant Q, the energy 
density decreases. As our goal is to assess the time at which the pressure obeys 
an equation of state (thereby justifying a hydrodynamical description of the 
system), the comparison of the relaxation for various couplings should be done 
for systems that have the same energy density. Therefore, in the comparison 
shown in fig. 11, the value of Q has been adjusted in each simulation such that 
the energy density remains unchanged. 

Fig. 11 demonstrates that the relaxation time decreases with increasing cou- 
pling constant g. In fig. 12 we have represented the relaxation time, defined 
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Figure 12: Points: relaxation time (see text for the definition used here) as a 
function of the coupling g. Line: fit by a power law. 



here as the time necessary to reduce the initial oscillations of the pressure by a 
factor 4, as a function of the coupling constant g for our set of values of g. One 
can fit all the points except the last one {g = 8) by a power law that suggests 
the following dependence^^ 

const 

The right most point in this plot is an outlier that does not follow this power 
law, possibly because this value of the coupling is too extreme for our approxi- 
mations/resummations to make sense (for g = 8, the interaction strength 5^/4! 
is significantly above 1). 



4.4 Energy density fluctuations 

The results we have shown thus far indicate that the pressure in the system 
relaxes towards the equation of state p — e/3, at relaxation times that decrease 
as the coupling constant increases. However, this study does not in and of itself 
tell us much about the nature of the state reached by the system. In particular, it 
does not tell us whether the system reaches a state of local thermal equilibrium. 
Because we have a system of strong fields whose modes have large occupation 
numbers, it is unlikely that the system can be described in terms of quasi- 
particles that have a Bose-Einstein distribution. In section 3, we observed that 

^^The axis of the figure 11 arc in lattice units. Thus, the horizontal axis is t/a and the 
vertical axis pa^ or where a is the lattice spacing. Since our model is scale invariant, 

the relaxation time scales like e"^/*. By eliminating a between the horizontal and vertical 
axis, it is easy to get the value of e^^^t. For g = 4 we have ea^ = 200 and the relaxation 
time is t/a 30, leading to e^/^t 113 (this combination is 11 for g = 8). Then, from one's 
favorite value of e in GeV/fm^, it is easy to obtain the relaxation time in fm's. 
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in the simple example studied that the phase-space density reaches a stationary 
form reminiscent of a micro-canonical equilibrium ensemble. Unfortunately, now 
that we are looking at a full fledged quantum field theory, the phase-space is 
infinite dimensional and whether the same behavior occurs is difficult to assess 
numerically. 

There are however signs of thermalization in the fluctuations of the energy 
distribution in the system. For the system as a whole, energy is conserved and 
will not fluctuate, regardless of whether the system is in thermal equilibrium or 
not. However, as is well known for canonical ensembles, by looking at energy 
fluctuations in a small subsystem, one can learn something about the energy 
exchanges between this subsystem and the rest of the system which acts as a 
heat bath. In particular, the nature of the fluctuations in the energy distribution 
of the subsystem can tell us whether it is in equilibrium with the rest of the 
system. If this is the case, the fluctuations are those of a canonical ensemble 
with a density operator p = e:iqi{—j3H). 

We show in figs. 13 and 14 results from a study for the smallest subsystem one 
can conceive of on a lattice-a single lattice site. In fig. 13, we display histograms 




Figure 13: Distribution of energy density at one lattice site, at various times in 
the evolution. The coupling constant is g — 0.5. 

of the values of the energy on one site^^ , at various times in the evolution. These 
curves are normalized so that their integral is unity-they can be interpreted as 
probability distributions for the value of the energy on one lattice site. At 
t = 0, this distribution is very close to a Gaussian, centered on the mean energy 
density in the system. The width of this Gaussian is entirely determined by the 
Gaussian spectrum of fiuctuations in eq. (25). At early times, the distribution 
first remains Gaussian-like, but tends to broaden with time. Around t « 30 in 
lattice units, we observe a rapid change of shape of this distribution-the peak 
of the distribution shifts to lower values of the energy and the tail extends much 

^''in lattice units, this is simply the value of T'^" at one given site. 
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further at large energy. Once this dramatic change of shape has taken place, 
the evolution of the distribution is rather slow and a stationary distribution is 
reached at late times. 

The evolution in the energy distribution can be explored further by looking 
at its moments defined by 

^ ■ (54) 

Higher moments are very sensitive to changes in the shape of the distribution, 
especially the appearance of an extended tail that signals broader energy fluc- 
tuations. We represent these moments as a function of time in fig. 14, up to 




Figure 14: Normalized moments of the energy density distribution 

at one lattice site, as a function of time. The coupling constant is g = 0.5. 

n = 6. They all start very close to 1 at i = 0, which is the sign of a very 
narrow distribution with little fluctuations. The rapid change of shape of the 
distribution around i « 30 corresponds to a rapid increase of the moments. By 
t ~ 70, the moments have reached nearly asymptotic values modulo moderate 
residual oscillations. 

It is interesting to compare the evolution of the energy distribution at a single 
lattice site with the time evolution of the pressure in fig. 9. The initial rapid 
decrease of the pressure oscillations is concomitant with the change of shape 
of the energy distribution. The subsequent (slower) relaxation of the residual 
oscillations of the pressure occurs after the energy has reached a stationary 
distribution. 



5 Summary and Outlook 

We discussed in this paper a formalism which resums secular terms in a weak 
coupling expansion of a scalar field theory with initial conditions generated by 
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strong sources. Wc showed that rcsummcd expressions, to aU orders in pertur- 
bation theory, for inclusive quantities could be expressed as an ensemble average 
of the corresponding leading order classical quantities where the initial classi- 
cal field for (;;idi member of the ensemble is shifted by a quantum fluctuation 
drawn from a Gaussian distribution. We showed that this averaging caused the 
resummed pressure to relax to a single valued relation with the energy den- 
sity and interpreted this as arising from the phase decoherence of individual 
classical trajectories. We showed in a toy model that for an expanding sys- 
tem our result leads to ideal hydrodynamical flow. We briefly addressed the 
issue of thcrmalization while our numerical results display features similar to 
those of a canonical thermal ensemble, they differ slightly in the particulars. A 
more systematic numerical study will likely be able to shed further light on this 
important point. As noted in the introduction, our system appears to satisfy 
Berry's conjecture which has been argued to be an important requirement for 
the thermalization of quantum systems. We plan to pursue this topic further in 
the future. 

Finally, we note that to be fully relevant to heavy ion collisions, our methods 
should be extended to gauge theories. We have shown previously that the 
formalism outlined in section 2 here is also applicable to a gauge theory [50]. 
The spectrum of quantum fluctuations [51] is the essential ingredient here and 
work on computing this quantity is well underway [52] . 

Acknowledgements 

We would like to acknowledge useful discussions with Jean-Paul Blaizot, Kenji 

Fukushima, Miklos Gyulassy, Tuomas Lappi, Larry McLcrran, Rob Pisarski, 
Andreas Schafer and Giorgio Torrieri. K.D's and R.V's research was supported 
by DOE Contract No. DE-AC02-98CH10886. F.G's work is supported in part 
by Agcncc Nationalc dc la Recherche via the programme ANR-06-BLAN-0285- 
01. We thank the Institute for Nuclear Theory at the University of Washington 
for its hospitality. One of us (FG) would like to thank Brookhaven National 
Laboratory as well as the Yukawa International Program for Quark-Hadron 
Sciences at Yukawa Institute for Theoretical Physics (Kyoto University) for 
partial support during the completion of this work. 

A Lattice implementation 

In our numerical calculation, we discretize space in a cubic lattice, while 
retaining time as a continuous variable. This means that when solving the 
classical equation of motion for the fleld, the timestep is freely adjustable in 
order to warranty a given accuracy. In this appendix, we summarize the main 
aspects of this lattice formulation. 
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A.l Discretization of space 

The field ip{t,x) (and its time derivatives) become a function of a continuous 
time t and of discrete indices i,j,k that vary between and L — 1. We im- 
pose periodic spatial boundary conditions to the fields. In order to keep the 
code simple, all the dimensionful quantities are expressed in lattice units, which 
amounts to choosing a lattice spacing a = 1. Thus, the classical equation of 
motion for the field becomes: 

'fijk{t) = iPi+ljk + fi-ljk + 'fiij + lk + fij-lk + ^^ijk+l + 'fiijk-l ~ Qfijk 

-Y^ljkit) + JiMt) ■ (55) 

(The terms on the right hand side of the first line correspond to the discretized 
version of the Laplacian of the field.) 

In our discretized description, there is an exactly conserved energy at > 0. 
To see that, multiply the previous equation by (fiijk and sum over all the lattice 
sites, 

di ^ + Ij^i^i+ljk-'Pijkf' + ■^{Vij+lk-^ijkf' + -^{Vijk+l-^ijkf 

ijk 

+ ^'Pijk\ = Yl Jijk'Pijk ■ (56) 

ijk 

The left hand side of this equation, which is nothing but the time derivative 
of the total lattice energy, is zero when the source J is turned off. But note 
that for this to work, we had to use a non-symmetric form of the discrete 
derivative in the definition of the kinetic energy. Taking the seemingly more 
natural ^(f i+ijk—'fii-ijk) instead oi f i+ijk—'fiijk would lead to an energy which 
is not exactly conserved on the lattice (though violations of energy conservation 
would be of higher order in the lattice spacing, it is preferable to use a lattice 
definition that makes it exactly conserved for any lattice spacing). 



A. 2 Small field fluctuations 

To obtain the spectrum of fluctuations for the initial condition of the classical 
field, we must also consider the time evolution from t = —00 to t = oi small 
perturbations to the classical field that behave like plane waves in the remote 
past. On the lattice, free plane waves are labelled by three integers I, m, n and 
are of the form: 

where the energy of the mode Imn is given by 

'27:l\ /2nm\ /27rn\\li/2 
— rns — rns 



2 I 3 — cos I j ~ cos y j ~ '^'-'^ y 



(58) 
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These plane waves obey the lattice Klein-Gordon equation, 

aijk = di+ljk + di-ljk + dij+lk + dij-lk + dijk+l + dijk-l — Sajjfe . (59) 

The lattice version of the fluctuations a±k{x) that enter in eqs. (25) are fluctu- 
ations a^j^'™"^(i) that obey 

.. 9^ 2 

^ijk + '^fijk ^ijk — 

= (Xj+lj/s + a,i-ljk + O-ii+lk + O'ij-lk + O-ijk+l + O-ijk-l ~ ^O'ijk ^60) 

and behave like the free waves of eq. (57) when t — )• — oo. 

A. 3 Sampling of the Gaussian fluctuations 

Once the fluctuations that obey this equation and initial conditions are known, 
the discrete version of eq. (25) reads 

mn 

and similar formulas for the correlators involving the time derivatives. Gen- 
crating Gaussian fluctiiations that have a given correlation function in general 
requires to diagonalize the 2-point correlator. On the lattice, this means that 
one should diagonalize an x matrix, a fairly time consuming step even for 
a reasonably sized lattice. 

In the case of the correlation function (61), it turns out that we can avoid 
this diagonalization. First of all, let us decompose a fluctuation a^fc on the 

basis formed by the a^jft , 

1 \ ^ 1 r (+lmn) r, { — Im-ajX I cin\ 

"'ij'' = 73 ^ijk + Plrnn « -jfe • (62) 

Imn 

Since we want to obtain a real fleld, we must have /3imn = C(*mn- Since there 
is a linear relation between aijk and the coefflcients aimn, l^imn, it is obvious 
that these coeflicients should also be Gaussian distributed. A little guesswork 
indicates that in order to obtain eq. (61), one needs 

(aimnai'm'n') — (PlmnPl'm'n') = , 

mm' ^nn' ■ (63) 

Equivalently, this can be translated into correlators for the real and imaginary 
parts of aimn- 

(Re {aimn)'Re {ai'm'n')) = (lm(aim„)Im («;/„/„/)) = ^^^^ Su'Smm'Snn' 

(Re (a;TO„)Im {ai'm'n')) = . (64) 
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In words, the real and imaginary parts of the coefScient aimn arc independent 
Gaussian random variables. This leads to a straightforward algorithm for gen- 
erating the correct Gaussian fluctuations of the initial conditions of the classical 
field: 

i. From t = — oo (in practice, some large negative time) to t = 0, solve the 
classical equation of motion (55) for the classical field cp given a source J, 

ii. For each Fourier mode lm,n, solve the equation of motion (60) for the small 

fluctuation until t = 0, 

iii. For each mode Imn, generate two random numbers Re (ajmn); Im {aimn) 
drawn from a Gaussian distribution of variance EimnL^/^, 



iv. Construct a field fluctuation a, a as 

(^ijk = 

^ijk = 7^ TT— Re aj„„ a|+j!'""\0) , (65) 



Imn 



aimnalj^ '{0} 



and superimpose it on to the fleld v'jjfc(O)) V'ufc(O) in order to obtain a 
fluctuating initial condition for the evolution at i > 0. Repeat steps iii 
and iv for each configuration in the Monte-Carlo evaluation of eq. (26) in 
order to obtain an ensemble of initial conditions for the classical field. 

V. For each initial condition constructed in this way, solve the classical equa- 
tion of motion (55) for t > (but now with J = 0). Compute the energy- 
momentum tensor of this classical field. Average it over the ensemble of 
initial conditions. 



A. 4 Ultraviolet sector 

The summation over fluctuations of the initial conditions for the classical field 
in eq. (26), if performed without any constraint on the momentum of the fluc- 
tuations one includes, leads to ultraviolet divergences in the energy-momentum 
tensor. 

For instance, if we impose a cutoff A on the momentum of the fluctuations 
included in eq. (26), one can check that generically the energy density will 
contain terms of the form 

e = ^ e g'A2 e a* , (66) 

.9 

where Q is the physical scale introduced via the source J. Eq. (26) is not 
renormalizable in the visual sense, because it results from a resummation that 
mixes diagrams with arbitrarily high loop orders, and it should be seen as an 
effective formula that only makes sense with an ultraviolet cutoff. In order to 



35 



include all the Fourier modes that are subject to instabilities, we must choose 
the cutoff such that Q < A. On the other hand, the cutoff should not be too 
large, otherwise the result of eq. (26) will be cutoff dependent. It turns out 
that at weak coupling {g ^ 1), there is ample room to choose such a cutoff 
(large enough to encompass the relevant physics and small enough to keep the 
cutoff-dependent terms small). Indeed, it is sufficient to have 

< A < . (67) 

This window of allowed A's can be enlarged if we notice that the A^ terms are 
pure vacuum contributions, that can computed and subtracted once for alP^. 
After this subtraction has been performed, the condition on A becomes 

g < A < ^ . (68) 
9 

Thus, taking A of order Q is a valid choice at small coupling. 



A. 5 Choice of the lattice cutoff 

We have seen in the previous subsection that a UV cutoff is necessary in order 
to ensure the finiteness of eq. (26). One should choose the cutoff so that all the 
modes up to the resonance band are comprised below the cutoff. We have also 
seen that if the coupling constant is small, then the dependence on such a cutoff 
is negligible since it occurs only in terms that are subleading by one power of 
9'. 

Such a cutoff exists naturally on the lattice, as a consequence of the dis- 
cretization. In the figure 15, we have represented the density of lattice Fourier 
modes as a function of their energy Eimn- This density falls abruptly when the 
energy reaches its maximal allowed value E = \fvl (in lattice units), and no 
Fourier mode with a larger energy can exist on the lattice. 

In lattice units (a = 1), the interplay between the physical scales and the 
lattice cutoff can be tuned by adjusting the parameter Q that sets the magnitude 
of the source J. The modes that are the most important for the decoherence 
responsible for the relaxation towards the equation of state arc the modes k < Q. 
If we chose a too large value of Q, the lattice cutoff will suppress modes that 
are important for this relaxation process. On the other hand, if Q is too low, 
then the physically rcilcvant modes fall in the region where the density of lattice 
modes (see the figure 15) is very small. In this case, very few lattice modes are 
available to represent the relevant physics. The optimal choice of Q is to take it 
so that the resonance band is located just before the fall off of the mode density 
at large Eimn- In this way, the physically relevant modes sit in the region where 
the lattice mode density is the largest. 

^''This is done by running the same simulation with the source J turned off, and by subtract- 
ing the corresponding result. We find that on a 12^ lattice, the vacuum contribution to T"" 
is almost independent of the coupling g, and equal to T^^c 1-35. The vacuum contribution 
to the pressure tensor is proportional to the identity, with a pressure equal to ^T°°^. 
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Figure 15: Density of Fourier modes as a function of their lattice energy (in the 
limit L — >■ oo). 



B Linear stability analysis for 0^ field theory 
B.l Instabilities of small perturbations 

The assumed suppression of the NLO correction relies on the fact that the 
perturbations /3,a±k introduced in section 2.3 have their "natural" order of 
magnitude {a±k ~ 0(l),/3 ~ 0{gQ)). If their equation of motion suffers from 
instabilities that amplify them, then the previous estimate is incorrect and the 
NLO corrections may be as large as the leading order contribution. It is therefore 
necessary to study the stability of small perturbations to the classical field (p. 

To keep the discussion simple, let us assume in this appendix that the source 
J, and therefore also the classical field do not have any spatial dependence^^. 
If a(x^ , x) is a perturbation to the field we can write its evolution equation 
in a linearized form as long as it remains small compared to the classical field 
(/? itself, 

This equation can be simplified by Fourier transforming the field a{x^ , x) in the 
spatial variables^^, 

a{x\x)^ j ^a{x',k)e^^-^ , (70) 

so that 

S+(fc' + yv')a = 0. (71) 

^^The results we obtain here remain valid in the case their spatial gradients are small. 
^^We use the same symbol for a and for its Fourier transform to keep the notations light, 
as the context always allows one to distinguish the two. 
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(The dot denotes a derivative with respeet to time.) Given a pair of solutions 
ai and of this equation, the Wronskian W = dia2 — aid2 is independent of 
time. 

When Lp depends only on time, it is a periodic function of time at .t" > 0. 
Therefore, the coefficient of the term a{x^ ,k) in eq. (71) is also periodic in 
time, which may lead to parametric resonance phenomena. For linear equations 
with periodic coefficients, the stability analysis can be performed by finding the 
"mapping at a period", that evolves a pair of solutions ai,2 from a;" = to 

=T (where T is the period of the coefficients in the equation): 

ai(r,fc) a2(r,fc)\ /ai(0,fc) 02(0, fc)\ 

di(T,fe) d2(T,fe)j Vdi(0,fe) d2(0,fe)j ■ ^^^^ 

This mapping can be written as a multiplication by a matrix Mfe thanks to the 
fact that equation (71) is linear. If the mapping Mfe is known, then after n 
periods one has 

ai(nT,fe) a2(nT,fc)\ /ai(0,fc) 02(0, fe)\ 

di(nT,fe) a2{nT,k))-^\ai{(),k) a2{0,k)J' 

and it is clear that the asymptotic behavior of the solutions ai^2 is determined 
by the eigenvalues Ai,2 of the matrix Mf,- 

From the conservation of the Wronskian, it is immediate to get 

det(Mfe) = A1A2 = 1 . (74) 

The two eigenvalues are thus mutually inverse, and we can write the trace as: 

tr(Mfe) = A + A-\ (75) 

where A is any of the two eigenvalues. One has therefore several cases: 

• If tr(Mfc) > 2 (resp. tr(Mfc) < —2), then A is real and greater than 1 
(resp. smaller than -1). In this case, the solutions of eq. (71) generically 
diverge exponentially with time. 



If tr(Mfc) ~ 2 (resp. —2), then A = 1 (resp. -1). The two eigenvalues 
of Mfc are in fact equal to 1 (resp. -1). For A = 1, this implies that the 
matrix is of the form 

In this case, one of the solutions is T-periodic (and is therefore stable), 
while the other solution of the basis diverges linearly with time (unless 
a = accidentally). 

If —2 < tr (Mfc) < 2, the eigenvalues Ai_2 are complex and lie on the unit 
circle, Ai,2 = exp(±z6') (they are mutual complex conjugates since the 
matrix is real valued). In this case, the small fluctuations are stable. 
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Figure 16: Trace of the monodromy matrix as a function of kmo (we denote 



One can compute numerically the monodromy matrix and its trace as a 
function of k. We have displayed the result in the figure 16. We see that 
the trace is equal to 2 for a discrete set of modes, including the zero mode^''. 
There is a band of modes^* 0.71 < k/mo < 0.76 in which the trace is greater 
than 2, and where the fluctuations are exponentially unstable. Outside of these 
discrete modes and of the resonance band, the trace is between -2 and 2, and 
the fluctuations are stable. Besides these rigorous statements about the location 
of the unstable modes, in practice the modes k for which |tr (Mfe)| < 2 but is 
close to 2 display a linearly growing behavior for a rather long time, before they 
eventually decrease. For instance, modes near the zero mode grow linearly for a 
time of order 27r/A; before they start decreasing. Strictly speaking, these modes 
are not unstable, but their growth at early times makes them important for the 
transient dynamics of the system. 



B.2 Lyapunov exponent 

In the resonance band, the exponential instability of the solutions of eq. (71) 
can be characterized by the Lyapunov exponent, that one can define from the 
largest eigenvalue of Mk as follows, 

^(fc. Too) = ^ In Max { Ai^s} , (77) 

^'^For fc = 0, it is easy to check that ai(t) = ip{t) and a2(t) = ip{t) dT /ifP{T) are solutions 
of eq. (71). Given these two solutions, it is straightforward to get tr (Mq) = 2. 

^°In section B.2, we show that the boundaries of this band are in fact I/V2 and 1/31/". 
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where jtiq = ^g^ipQ. Asymptotically, the solutions a{x^,k) of eq. (71) grow like 

The Lyapunov exponent can be obtained numerically from the trace of Mfc, but 
in the case of the (f)^ potential it can in fact be derived analytically^^. Consider 
the equation 

a + (fe^+m2(f))a = 0, (79) 

where 

1 

m + -TO^ = . (80) 
o 

(In other words, m{t) = g(p{t) / \/2- The equation for m{t) is a consequence of the 
equation for ip{t).)^^ Let us call mo the maximal amplitude of the oscillations 
of m{t), and introduce the new variable z defined by 

TO^ = mlz . (82) 

The time t and the variable z are related by 



f =mo^^.(l-.^), (83) 

and we can rewrite eq. (79) as 

2z{l - z'^)a" + (1 - 3z^)a' + 3{k^ + z)a = 0, (84) 

where the prime denotes a dcirivative with respect to z and k = k/mo. By this 
transformation, we have turned an equation with oscillating coefficients into an 
equation with polynomial coefficients. Given a pair ai_2 of solutions of eq. (84), 
one can show that the Wronskian is 

W = a[a2 - aia'2 = ^" , (85) 
^/z{l - z^) 

where wq is a constant. 

Let us call now M = aia2, with ai,2 two solutions of eq. (84) (possibly 
identical). A straightforward calculation shows that M obeys the following 
third order differential equation: 

2z{l - z^)M"' + 3(1 - 3z'^)M" + 6{z + 2k^)M' + 6M = . (86) 



^'^Thc derivation we expose here is a particular case of techniques developed in [30]. 
^''The value of g is not important per se. Only the combination mo ~ gipo matters for the 
Lyapunov exponent. Moreover, if ;u(fc,mo) is the Lyapunov exponent, then it scales as 

VA , /u(Afe, Amo) = A;u(fe, mo) , (81) 

due to the scale invariance of our model. 
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If ai and 02 arc two independent solutions of eq. (84), then the three independent 
solutions of eq. (86) can be thought of as , a| and aia2- A remarkable property 
of eq. (86) is that it admits a polynomial solution: 



M{z) 



2k^z + Ak^ ■ 



1 . 



(87) 



If K is in the resonance band, where ai increases exponentially while 02 de- 
creases exponentially, this polynomial solution must be their product M{z) = 
ai{z)a2{z). With the help of the Wronskian, it is then easy to find 



01(2;) = y/\M{z)\ exp 



' 2 



dz 



a2{z) 



exp 



M{z)^z{l - z2) 
dz 



M{z)^Jz{l- z"^) 



(88) 



In order to determine the constant wq, one must insert these solutions into 
eq. (84). This leads to 



Wo = 4W6k2 



(89) 



The resonance band corresponds to the values of k such that the argument of the 
square root is positive (otherwise wq would be imaginary and one would have 
an oscillating solution instead of an exponentially growing solution). Thus, the 
instability domain is 

^ < K < -TTj . (90) 
^/2 ~ ~ 31/4 

From the above solutions ai, 02, one can also determine their growth during one 
period of oscillation of m(t), from which one gets the Lyapunov exponent. One 
obtains 

dz 



IJ,{k,mo)T = 2wo / 

Jo 



M(z)^z{l~z^) ' 
where T is the period of the oscillations of m(t): 

4V6 



T = 



mo 



We finally get: 



IJ,{k, mo) = 2mo4 / 



3 " 



Jo 7T^' 

ri dz 

l\JO (z2_2k2z+4k4_1)^^(1_22) 



Jo 



(91) 



(92) 



(93) 



The integrals in this formula can be evaluated numerically^^ , and one can com- 
pare this analytical result to the direct numerical computation (performed by 
computing the trace of Mk). As illustrated in the figure 17, the agreement is 
perfect. 



^^The integral in the numerator must be handled as a Cauchy principal value, since its 
integrand has two poles in the interval z 6 [0,1]. 
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Figure 17: Comparison between a numerical computation of the Lyapunov ex- 
ponent, and its evaluation from eq. (93). 

B.3 Linear instability and time decoherence 

Note that a zero Lyapunov exponent does not mean that the solutions are stable, 
it only implies that they are not exponentially unstable: outside of the resonance 
band determined by the inequality (90), the solutions of eq. (71) may exhibit a 
linear growth in time - either indefinitely if tr (M^) = 2 or during a finite time 
for all the other modes. In the case of the zero mode, this linear behavior has 
a very simple interpretation, because it is a direct consequence of the fact that 
the (j)'^ potential leads to non-harmonic oscillations. Indeed, the classical field (p 
oscillates periodically in the potential with a frequency uj that is proportional 
to the amplitude ipo of the oscillations'^^: a classical field that oscillates freely 
in a (j)'^ potential can be written as 

^{t) = <fiof{^ot) , (94) 

where / is a periodic function of amplitude unity. Let us now add a small 
perturbation a to this classical field, so that its amplitude is now ipo + uq. The 
perturbed oscillations are given by: 

^it)^{ipo + aQ)fii^o + ao)t) . (95) 

If we assume that oq ^ (po and we linearize in the perturbation, we have 

i;{t) = p{t) + a^fipot) -t- (^0 aot f{pot) + 0{al) . (96) 

^^This result is specific to a potential that has only a rf>^ term. Indeed, in four dimensions, 
such a field theory is scale invariant at the classical level, and its only dimensionful parameter 
is the amplitude of the oscillations of (p (set by the external source J that drives the system 
at < 0). For other potentials, the oscillation frequency lj is in general some complicated 
function of the amplitude ipo and of the coupling constants present in the potential. 
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The fact that the frequency depends on the amphtude produced a term that 
grows hnearly with time for the hnearized perturbation. This result is ubiquitous 
for any non- harmonic potential, and not specific to (f)^. This is the reason why 
a linear instability is the generic behavior of the solutions of eq. (71). Due to 
its origin, this linear instability is also closely related to another property: two 
classical solutions that at = differ only by a small perturbation ao will 
eventually become completely incoherent since their phase has shifted by 2tt 
after a time of order 

^decoherence ^ • ) 

In the case of perturbations that have a non- zero momentum k, the scaling law 
(94) is not valid if g<fo < |fe|. For these high fe modes, the oscillations are almost 
harmonic and thus independent of their amplitude. This means that the linear 
growth becomes weaker and weaker as |fe| increases, and practically irrelevant 
for Fourier modes above gipo. 
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